clear all
*cap log close
set more off
set seed 603

* show IV estimate hetereogenity by motorist characteristics ----------------

* setup data with characteristics ------------
use "${data}/out/4-main", clear
gen Y = cite_ny1 

* heterogeneity groups -----------
gen old = age>30
label var old "Age > 30"
gen white = (race=="W")
replace white = . if race=="U"
label var white "Race = White"
gen rich = zipincome > 50
label var rich "Zip Income > 50K"
replace rich = . if zipincome_miss==1
gen pastyear = cite_py1
label var pastyear "Citation Past Year"
* ----------------------------------------------


* compute + store IV estimates -----------------
gen ind = .
gen lab = ""
gen val = .
gen m   = .
gen b   = .
gen s   = . 

* loop over het dimensions --------
local k = 1
foreach x in pastyear rich local white old female {
	replace ind = `k'   if _n == `k'
	replace lab = "`x'" if _n == `k'
	replace val = 0     if _n == `k'
	qui summ Y if `x'==0 & lenient == 1
	replace m = `r(mean)' if _n == `k'

	ivreghdfe Y (harsh = Z) if `x'==0, absorb(totfe) vce(cluster officerid)
	replace b =  _b[harsh] if _n == `k'
	replace s = _se[harsh] if _n == `k'
	local ++k 

	replace ind = `k'-1 if _n == `k'
	replace lab = "`x'" if _n == `k'
	replace val = 1     if _n == `k'
	qui summ Y if `x'==1 & lenient == 1
	replace m = `r(mean)' if _n == `k'

	ivreghdfe Y (harsh = Z) if `x'==1, absorb(totfe) vce(cluster officerid)
	replace b =  _b[harsh] if _n == `k'
	replace s = _se[harsh] if _n == `k'
	local ++k
}

* get treatment effects relative to the mean ---------
gen te = m + b 
gen ub = m + b + 1.96*s
gen lb = m + b - 1.96*s

* setup for spacing in plot ---------------------
gen newind = ind 
replace newind = ind + 0.5 if val == 1
// -------------------------------------------------------

* store overall control mean for plot --------
qui summ Y if lenient == 1
local mu_ovr = `r(mean)'



// construct plot ---------------------------------------
#delimit ;	
twoway 
	(rbar ub lb ind if val ==0 , yaxis(1 2) horizontal barwidth(.3) color(pink%10)) // bar is CI: group1
	(scatter ind m if val ==0,  yaxis(1 2) mcolor(pink) msymbol(Oh) msize(large)) // dot is mean: group1
	(pcarrow ind m newind te if val ==0,  yaxis(1 2) lcolor(pink) mcolor(pink) lwidth(medthick)) // arrow end-y end-x arrowhead-y arrowhead-x: group1
	(pcarrowi 11 0.399 11 .383 (3) "Lenient Mean", lcolor(gray) mcolor(gray) mlabcolor(gray) barbsize(tiny)) //arrow end-y end-x arrowhead-y arrowhead-x (clock position of text)
	
	(rbar ub lb newind if val ==1 , yaxis(1 2)  horizontal barwidth(.3) color(blue%10)) // bar is CI: group2
	(scatter newind m if val ==1, yaxis(1 2) mcolor(blue) msymbol(X) msize(large)) // dot is mean: group2
	(pcarrow newind m newind te if val ==1, yaxis(1 2) lcolor(blue) mcolor(blue) lwidth(medthick)) // arrow end-y end-x arrowhead-y arrowhead-x: group2
	(pcarrowi 8.1 .425 8.8 .415 (5) "IV Estimate", lcolor(gray) mcolor(gray) mlabcolor(gray) barbsize(tiny)) //arrow end-y end-x arrowhead-y arrowhead-x (clock position of text) 

	, legend(row(1) order(6 "= 1" 2 "= 0") pos(6) region(lstyle(border))) //group markers' legend
	ylabel(1.25 "Citation Past Year" 3.25 "Zip Income > 50K" 5.25 "County Resident" 7.25 "Race = White" 9.25 "Age > 30" 11.25 "Female", angle(0) noticks nogrid axis(1)) //label yaxis left
	ylabel(1.25 "Citation Past Year" 3.25 "Zip Income > 50K" 5.25 "County Resident" 7.25 "Race = White" 9.25 "Age > 30" 11.25 "Female", angle(0) noticks nogrid axis(2)) //label yaxis right
	ytitle("") //suppress y title
	xline(`mu_ovr', lpattern(dash) lcolor(gs8)) // add dotted line
	xtitle("")
	xlabel(.2(.1).5,nogrid)
	graphregion(color(white)) bgcolor(white) plotregion(lcolor(black) lwidth(medthin)) ;
#delimit cr
graph export "${out}/apx_deter/heterogeneity.pdf", replace




